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Abstract 


A time- accurate approximate- factorization (AF) algorithm is described for solution of the three- 
dimensional unsteady transonic small-disturbance equation. The AF algorithm consists of a time- 
linearization procedure coupled with a subiteration technique. The algorithm is the basis for the 
CAP-TSD (Computational Aeroelasticity Program-Transonic Small Disturbance) computer code, 
which was developed for the analysis of unsteady aerodynamics and aeroelasticity of realistic aircraft 
configurations. The paper describes details on the governing flow equations and boundary conditions, 
with an emphasis on documenting the finite-difference formulas of the AF algorithm. 

Introduction 

Considerable research is being conducted presently to develop computational fluid dynamics 
methods for predicting unsteady transonic aerodynamics for aeroelastic applications (ref. 1.) The 
resulting computer codes are being developed to provide accurate methods of calculating unsteady 
air loads for the prediction of aeroelastic phenomena such as flutter and divergence. One of the 
most fully developed codes for analysis of transonic aeroelasticity, for example, is the CAP-TSD 
(Computational Aeroelasticity Program-Transonic Small Disturbance) computer code (ref. 2.) The 
code permits the calculation of unsteady flows about realistic aircraft configurations for analysis of 
aeroelasticity in the flutter-critical transonic speed range. It can treat configurations with general 
combinations of lifting surfaces and bodies. 

The code uses a time-accurate approximate-factorization (AF) algorithm (refs. 3 and 4) developed 
for solution of the three-dimensional unsteady transonic small-disturbance (TSD) equation. The AF 
algorithm involves a time-linearization procedure coupled with a subiteration technique and is similar 
to the unsteady full-potential algorithm reported by Shankar et al. (ref. 5.) References 3 and 4 show 
the AF TSD algorithm is efficient for application to steady or unsteady transonic flow problems. It 
can provide accurate solutions in only several hundred time steps to yield a significant computational 
cost savings compared with alternative methods. This paper describes details of the governing flow 
equations and boundary conditions, with an emphasis on documenting the finite-difference formulas 
of the AF algorithm that were not reported in references 3 and 4. 

Symbols 

A coefficient in TSD equation (see eqs. (3)) 

B coefficient in TSD equation (see eqs. (3)) 

C coefficient in nonreflecting far-field boundary conditions (see eq. (14a)) 

C a angle-of-attack correction for bodies (see eqs. (12a) and (12b)) 

C p pressure coefficient 

Cp, isentropic pressure coefficient 

C Ps pressure coefficient due to change in entropy 

C Pv pressure coefficient due to vorticity 

Ct thickness correction for bodies (see eqs. (12a) and (12b)) 

c v specific heat at constant volume 

D coefficient in nonreflecting far-field boundary conditions (see eq. (14b)) 

E coefficient in TSD equation (see eqs. (3)) 

F coefficient in TSD equation (see eqs. (4) to (7)) 

F\ flux in L £ operator (see eqs. (69) and (70)) 
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flux in Ly operator (see eq. (95)) 
flux in Lq operator (see eq. (Ill)) 
flux in time derivative in TSD equation (see eq. (2a)) 
flux in x-direction in TSD equation (see eq. (2b)) 
flux in y-direction in TSD equation (see eq. (2c)) 
flux in ^-direction in TSD equation (see eq. (2d)) 
coefficient in TSD equation (see eqs. (4) to (7)) 
flux in time derivative in TSD equation (see eq. (132a)) 
flux in ^-direction in TSD equation (see eq. (132b)) 
sonic reference flux (see eqs. (144) and (146)) 
flux in 77-direction in TSD equation (see eq. (132c)) 
flux in ^-direction in TSD equation (see eq. (132d)) 
coefficient in TSD equation (see eqs. (4) to (7)) 
index in £- or x-direction 

index in 77- or y-direction 

__ 1 

1 + Q 

index in or z-direction 
differential operator in ^-direction 
differential operator in 77-direction 
differential operator in ^-direction 
free-stream Mach number 
component of normal vector on body 
component of normal vector on body 
component of normal vector on body 
number of grid points in ^-direction 
number of grid points in 77-direction 
number of grid points in ^-direction 
index denoting time level 
constant (see eq. (16a)) 
residual (see eq. (131)) 

local change in entropy from free-stream value 
nondimensional time 
shock speed 

flow speed upstream of shock (see eq. (20)) 
speed term (see eq. (16c)) 


in x-direction (see eq. (11)) 
in y-direction (see eq. (11)) 
in ^-direction (see eq. (11)) 
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V s sonic reference speed (see eq. (16b)) 

W defined in equation (84) 

x physical coordinate in streamwise direction 

y physical coordinate in spanwise direction 

z physical coordinate in vertical direction 

at angle of attack of body 

/ % yaw angle of body 

T wake circulation 

7 ratio of specific heats 

At nondimensional time step 

Ati step size from time level n to time level n + 1 (see eq. (65)) 

At 2 step size from time level n- 1 to time level n (see eq. (66)) 

A*3 step size from time level n - 2 to time level n- 1 (see eq. (67)) 

C computational coordinate in vertical direction 

7] computational coordinate in spanwise direction 

£ computational coordinate in streamwise direction 

a relaxation parameter (see eq. (31)) 

potential function (see eq. (23)) 

(f> disturbance velocity potential 

(j) first intermediate disturbance velocity potential 

<f> second intermediate disturbance velocity potential 

4>x disturbance velocity in x-direction 

<i py disturbance velocity in y-direction 

<f>z disturbance velocity in 2-direction 

function representing stretching and rotating of vortex filaments associated with 
entropy variation (see eq. (23)) 

Governing Equations 

TSD Equation 

The flow is assumed to be governed by the general-frequency, modified, unsteady TSD potential 
equation, which may be written in conservation law form as 


where 


dfo,dh,df 2 ,dh =n 

dt dx dy dz 

(1) 

fo = -Afa - B<t> x 

(2a) 

h = E<t> x + F4> 2 x + G<t>l 

(2b) 

h — 4>y + H4> x 4>y 

(2c) 

II 

(2d) 
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The coefficients A, B , and E are defined as 


A = M% 0 B = 2M% 0 and E = \-M% ) (3) 

Several choices are available within CAP-TSD for the coefficients F, G, and H depending upon the 
assumptions used in deriving the TSD equation. Briefly, the coefficients are referred to as “NASA 
Ames” coefficients when defined as 

F = -i( 7 + l)M£ G=±( i-mlo H = -( n -l)M 2 00 (4) 

and are referred to as “NLR” coefficients when defined as 

F = ~\ [ 3 - ( 2 - 7)^] Ml G = -\mI H = -Ml (5) 

The “classic” coefficients are given by 

F = ~\{l + mL ^=0 tf=0 (6) 

and finally the coefficients for the linear equation are 

F = G = H = 0 (7) 


Coordinate Transformation 

The unsteady transonic small-disturbance equation is solved numerically on a finite-difference 
grid in a computational coordinate system (£, 77, (). The finite-difference grids in both the physical 
and computational domains are contained within rectangular boundaries and conform to the leading 
and trailing edges of the horizontal lifting surfaces. Regions in the physical domain are mapped into 
rectangular regions in the computational domain using shearing transformations. For simplicity, no 
shearing is performed in the vertical direction, so that pylons and vertical tails are approximated by 
rectangular surfaces. The shearing transformation may be written generally as 


£ = £(s.y) 1 = y C = 2 


( 8 ) 


where f , 77, and C are the nondimensional computational coordinates in the x-, y-, and ^-directions, 
respectively. The TSD equation (eq. (1)) may then be expressed in computational coordinates as 


~dt + + a| + + G (£s + <Pt ,) 1 + ^ + (j) V ) + + <Py) 

+ ^n h + ^ + + ^(^' ?;>c ) =0 


(9) 


Boundary Conditions 

The horizontal lifting surfaces are modeled (ref. 6) by imposing the following boundary conditions. 
Flow tangency: 


4t = g + ft 


(10a) 



Trailing wake: 


[ 0 2 ] = 0 

[<t>x + <f>t] = 0 


(10b) 

(10c) 


where / is a function describing the position of the lifting surface (including thickness, camber, angle 
of attack, yaw angle, and dihedral) and [ ] indicates the jump in the indicated quantity across the 
wake. The flow-tangency condition is imposed along the mean plane of the respective lifting surface. 
In equation (10a), the plus and minus superscripts indicate the upper and lower surfaces of the mean 
plane, respectively. The wakes are assumed to be flat and horizontal. 

Bodies such as the fuselage, stores, and nacelles are treated as follows. For a body at angle of 
attack a b and at yaw angle 0 b , the exact steady flow-tangency boundary condition may be written 
as 

N X ( 1 + 0x) + Ny((fiy + (3 b ) + N z ((pz + cq>) = 0 (11) 

where N(x , y, z, t) = 0 defines the body surface. Computationally, bodies are modeled with simplified 
boundary conditions applied on a computational box that extends to the upstream and downstream 
boundaries of the grid with a rectangular cross section rather than on the true surface (refs. 7 
and 8.) The method is consistent with the small-disturbance approximation and treats bodies with 
sufficient accuracy to obtain the correct global effect on the flow field without the use of special 
grids or complicated coordinate transformations. As such, the approximations to the flow-tangency 
boundary condition (eq. (11)) imposed on the computational box are 

0± = -C't(^) ± -C a /3 6 (12a) 

for right or left surfaces and 

4>t = ~ c t ~ Can, > ( 12b ) 

for top or bottom surfaces. The parameters Ct and C a are thickness and angle- of- attack corrections, 
respectively, derived from slender-body theory to account for spatial differences between true and 
computational body surfaces (refs. 7 and 8.) Also, in equation (12a), the plus and minus superscripts 
indicate the right and left surfaces, respectively. In equation (12b), the plus and minus superscripts 
indicate the top and bottom surfaces, respectively. 

The conditions imposed upon the outer boundary of the computational region are similar to the 
characteristic or “nonreflecting” boundary conditions reported by Whitlow (ref. 9.) The conditions 


employed here 

are given by the following. 


Upstream: 

0 = 0 

(13a) 

Downstream: 


(13b) 

Above: 

+ <t>z= 0 

(13c) 

Below: 

-<t>z = o 

(13d) 
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Right spanwise: 


D 

+ 4>y = 0 

(13e) 

Left spanwise (for full-span modeling): 


D 

— 4>y = 0 

(13f) 

Symmetry plane (for half-span modeling): 


(j)y = 0 

(13g) 

where 


C = E + 2 F(j) x 

(14a) 

/ 1/2 
D -[ 4 A+ o) 

(14b) 


Entropy Model 

An entropy model was developed for the AF algorithm (ref. 10), similar to that of reference 11, 
that includes (1) an alternative streamwise flux, (2) an entropy correction in the pressure formula, 
and (3) a modified wake boundary condition to account for convection of entropy. In this section, 
the entropy model is described. 

Alternative streamwise flux. The entropy model is formulated by first replacing the 
streamwise flux f\ (eq. (2b)) in the TSD equation with an alternative flux given by 

fl = (7 + 1 )MlQ (vv s - l -v 2 ) + G4>1 (15) 

where 


[2 + (7 — 1)11*2,1 1/2 
[ (7 + 1^ 

(16a) 

Y s - ^ " 1 
2 Q 

(16b) 

y _ (1 + Q)<f>X 

1 + <f>x + Q 

(16c) 

Pressure correction. The pressure formula is modified to include entropy effects according to 

C p — C Pi + Cp s 

(17) 

where C Pi is the isentropic pressure coefficient and C Ps is the pressure coefficient due to the change 
in entropy. As reported in reference 10, C Ps is given by 

c — 2 s 

Ps ~ 7(7 - Cv 

(18) 
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where s is the change in entropy from the free-stream value. Equation (18) obviously requires the 
determination of entropy along the surface of the geometry being considered. This first requires 
the determination of shock location and then the calculation of the entropy jump across the shock. 
The shock location is easily determined since most TSD algorithms use type-dependent switching 
to capture shocks and to properly treat regions of subsonic and supersonic flow. The entropy jump 
is computed with the Rankine-Hugoniot shock jump relation: 


s_ (7 + l)«i - (7 - i)Q 2 
Cy n (7 + l )Q 2 - (7- !) u i 


1 

- 11 "? 


(19) 


where 

u\ = 1 + (j) x - u s (20) 

In equation (20), u\ is the flow speed upstream of the shock and u s is the shock speed. In the present 
formulation, the entropy is convected downstream from the shock according to 


ds ds _ 
dt dx 


(21) 


The correction to the pressure formula to include entropy effects (eq. (17)) does not directly affect 
the flow field. The effect on the flow field is produced by the modified wake boundary condition. 


Modified wake boundary condition . The wake boundary condition requires that the pressure 
be continuous across the wake. Since the pressure formula (eq. (17)) includes a term due to entropy, 
the wake boundary condition must be similarly modified as 

r* + r x = - (22) 

where A represents the jump across the wake. In equation (22), A C Ps is determined by first 
convecting the entropy along the wake and then computing Cp s with equation (18). The nonzero 
right-hand side of equation (22) thus modifies the circulation distribution T. 

Vorticity Model 

A vorticity model was developed for the AF algorithm (ref. 10) similar to that of reference 12. 
In this section, the vorticity model is described in detail, including (1) a modified velocity vector 
that in turn modifies the TSD equation, (2) a pressure formula correction for vorticity effects, and 
(3) the resulting wake boundary conditions. 


Modified velocity vector . The vorticity model is formulated by first writing the velocity vector 
as the sum of potential and rotational parts according to 


U = V$ i—— v* 

7~ lc i> 


(23) 


In equation (23), the first term is the gradient of a scalar potential <E> and the second term is the 
product of the change in entropy s and the gradient of the function '!' . The function 4/ is a measure 
of the stretching and rotating of vortex filaments associated with entropy variation (ref. 12). In the 
present algorithm, the rotational part of the velocity vector is assumed to occur only in the region 
downstream of shocks. Further assuming that the entropy convects with the free stream and that 
the shock curvature is negligible implies that 


d'l' _ 1 _ o d$_ 

dx J^oo 


(24) 
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These assumptions eliminate the variable from the model and leave only the change in entropy to 
be determined throughout the flow field. In a steady flow, entropy is constant along streamlines and 
changes only through shock waves. The entropy jump is computed along shocks with the Rankine- 
Hugoniot relation (eq. (19)). Then, for simplicity, the grid lines are assumed to approximate the 
streamlines of the flow, which is consistent with the small-disturbance approximation. The entropy 
is then convected downstream along the grid lines for unsteady applications (eq. (21)) or is held 
constant along the grid lines for steady applications. 

The modified velocity vector in turn modifies the TSD equation because the streamwise distur- 
bance speed u = (f) x is now given by 


^ — 4>x 


1 s 

7(7 - 1 ) M$ Q Cv 


(25) 


The new TSD equation has the same conservation law form as equation (1), with new fluxes defined 
by simply replacing <f> x by the modified speed given in equation (25). 


Pressure correction. The pressure formula must also be modified when vorticity effects are 
included in the model. In general form, the pressure coefficient may be computed from 

Cp = C Pi + C Ps + C Pv (26) 

where C Pv is the pressure coefficient correction due to vorticity. The correction due to vorticity 
approximately cancels the correction due to entropy, and thus the pressure coefficient C p is given by 
the isentropic formula. At the TSD equation level, this is clearly demonstrated by first considering 
the general form of equation (17). Assuming the first-order small-disturbance pressure formula for 
C Pi , defining C Ps as given by equation (18), and replacing <f> x by the modified disturbance speed of 
equation (25) yields 


C p 2 '4t 2<p x 7(7 _ 1)M 2 3Cu + 7(7 _ 1)M 2 )Cu (27) 

Here the corrections due to entropy and vorticity identically cancel each other, and thus the pressure 
coefficient is given by the isentropic formula in terms of the irrotational disturbance speed <j> x . 

Modified wake boundary condition. As with the entropy model, the wake boundary condition 
in the vorticity model requires that the pressure be continuous across the wake. Since the pressure 
is now given by the isentropic formula (eq. (27)), the wake boundary condition is identical to the 
original condition given by 

r t + r x = o (28) 

Consequently, the modifications that are required when both entropy and vorticity effects are 
included are the alternative streamwise flux given by equation (15) and the modified streamwise 
disturbance speed given by equation (25). 

Approximate-Factorization Algorithm 

General Description 

The AF algorithm consists of a time-linearization procedure coupled with a subiteration tech- 
nique. For unsteady flow calculations, the solution procedure involves two steps. First, a time- 
linearization step (described below) is performed to determine an estimate of the potential field. 
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Second, subiterations are performed to minimize linearization and factorization errors. Specifically, 
the TSD equation is written in general form as 

R (</>” +1 ) = 0 (29) 

where 0 n+1 represents the unknown potential field at time level n + 1. The solution to equation (29) 
is then given by the linearization of equation (29) about 0*: 

w+ (t?L. A *=° <3o) 

In equation (30), 0* is the currently available value of 0 Tl+1 and A 0 — 0 n+1 — 0*. During convergence 
of the iteration procedure, A 0 approaches zero so that the solution is given by 0 n+1 = 0*. In general, 
only one or two iterations at a given time level are required to achieve acceptable convergence. For 
steady flow calculations, iterations are not used since time accuracy is not necessary when marching 
to steady state. 

Mathematical Formulation 

The AF algorithm is formulated by first approximating the time derivative terms (<j>tt and (j> x t) 
by second-order-accurate finite-difference formulas. The TSD equation is rewritten by substituting 
0 = 0*4* A 0 and neglecting squares of derivatives of A0, which is equivalent to applying equation (30) 
term by term. The resulting equation is then rearranged and the left-hand side is approximately 
factored into a triple product of operators yielding 

L^L C A 4> = -aR (f, 4> n , <t> n ~\4> n ~ 2 ) (31) 

where the L £, and L ^ operators and the residual R are defined and described in subsequent 
sections. In equation (31) a is a relaxation parameter that is normally set equal to 1.0. To accelerate 
convergence to steady state, the residual may be overrelaxed using a > 1. Equation (31) is solved 
with three sweeps through the grid by sequentially applying the operators as follows. 


£-sweep: 

A 0 = ~aR 

(32a) 

r/- sweep: 




A 0 = A 0 

(32b) 

C-sweep: 




A 0 = A 0 

(32c) 


Time-Linearization Step 

An initial estimate of the potentials at time level n4-l is required to start the subiteration process. 
This estimate is provided by performing a time-linearization calculation. The equations governing 
the time-linearization step are derived in a similar fashion as the equations for subiteration. The 
only difference is that the equations are formulated by linearizing about time level n rather than 
about the iteration level. This is accomplished by substituting 0 = 0 n 4- A 0 into the TSD equation 
(eq. (1)) and neglecting squares of derivatives of A 0, as done previously. 

Difference Equations for the Disturbance Velocity 

The AF algorithm is simplified greatly, both mathematically and numerically, by first determining 
the disturbance velocity components 0*, 0*, and 0*. These components are required in numerous 
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places throughout the AF solution procedure at the half-node points (i±V 2 , j , A;), (i, j±V 2 , A;), and 
(i, j, fc±V 2 ). The finite-difference formulas that approximate the velocity components are presented 
in this section. In these formulas the i, j, or A; subscripts are often omitted for clarity. Also, the grid 
metrics that appear in the formulas are defined in the following section. 


At Half-Node Points (z±V 2 , j, fc) 

All the expressions for the velocity components are straightforward central-difference approxima- 
tions that are centered at the half-node points to be consistent with the treatment of the fluxes in 
the AF algorithm. For example, at half-node points (i— ' V 2 , j, k), the x- and y-components of the 
disturbance velocity are required. They are defined simply by 


^ x i-l/2j ^ X i~ 


1 

1/2 ' j 6-fc-i 


(33) 


*Vi- 1/2J 


. lziLj. 4 + 

, yi-l/2J £ _ 2 - TJj-l) 


If vorticity effects are included, then the x-component is defined by 


(34) 


< ^ x i-l/2,j ^ x t- 


~ ff-i 

1/2d & - &-i 27(7 


AufT [(-)* + (“)* 1 (35) 

-i)Mx>IA<Wi Vcv/j-iJ 


Formulas for the components at (i+V 2 , j , fc) are determined from equations (33) to (35) by 
incrementing the i index by one. Also, in the symmetry plane (normally taken to be j = 1) 
the y-component of the disturbance velocity is set equal to zero to impose the symmetry condition 
(eq. (13g)). Namely, 

= 0 < 36 > 

Furthermore, the y-component of the disturbance velocity must be defined differently to account for 
the flow-tangency condition (eq. (12a)) for bodies. The side surface of the computational box that is 
used to model a body, for example, is located an equal distance between grid planes in the spanwise 
direction (i.e., between planes j and j — 1). For grid points in the plane j , which is adjacent to 
the side surface, the ^-component of the disturbance velocity is defined as the weighted average of 
values at j+V 2 (determined by finite differencing) and at j—V : 2 (given by the boundary condition). 
The resulting formula is 


,* _ Vj-Vj-i 1 / ; 


Kj± 1 ~ ^i-l J+l 


& “ 1 


+ ' 


& - 6-1 


+ 


<t>l 


«J+1 


Vj+1 - Vj 


z3a 1 Wi 


n- u' 


T lj + 1 Vj 


Vj+1 - Vj 
Vj + 1 - Vj - 1 




(37) 


Similarly, for a full-span configuration, the ^-component of the disturbance velocity is defined 
differently to treat the left side surface (located between planes j and j + 1). Hence, for grid 
points in the plane j that is adjacent to the left side surface, 
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( 38 ) 


<t> 


_ Vj+1 Vj 1 < h 2 j 


< Ki-\ - fl-ij-i 


^-V2,i - w - Vj . x 2 V 6-6-! + ^'- x 6~6-l 


, *hz3±± | Cu-Cu-l ' 

Vj - Vj - 1 Vj - Vj-1 


JUi 

Vj+1 - Vj 


'■^h h (t) ~ + *■* + c “-‘ (t)L + C “‘- A 


At Half-Node Points (i, j±V2, fc) 

At half-node points ( i , j— V2, fc) the x- and ^-components are defined by 


,* * fli+lj ~ ^i-lj + ^-fl,j-l ~ ffi-lj-l 

^.i-V2 2(6+i“6-l) 

. = , + i-Cu-i , 

%-l/2 ^Vij-1/2 2(6+1 ~6-l) Vj ~ Vj-1 

If vorticity effects are included, then the x-component is defined by 

- faj-Cu + Wi-Cu-i 
9x U~ W hl ~ 1 ' 2 2 (^+ 1 - 6 - 1 ) 


27(7-1 )M2, [Vcvjj + 


(39) 

(40) 


(41) 


Formulas for the components at (z, J+V2, A:) are determined from equations (39) to (41) by 
incrementing the j index by one. Also, the formulas need to be modified to again account for 
the symmetry condition (j = J; normally taken to be J = 1). At j = J, the x- and ^-components 
become 


1/2 “ ^k,j+ 1/2 

^Vi.J-l/2 ~ ~^Vi,J+ 1/2 


(42) 

(43) 


Furthermore, both the x- and ^-components of the disturbance velocity must be defined differently 
to account for the fiow-tangency condition (eq. (12a)) for bodies. For half-node points (z, j— V 2, k) 
that lie on the right side surface of the computational box, a one-sided formula is used for the 
x-component given by 


,* _ 1 / Vj ± 1 - Vj-1 c fil+hj ~ ti-lj _ Vj - T?j-1 ffs+lj+l ~ ^-l,j+l \ 

2 y Vj+i — Vj 6+1 -6-1 7 J+ 1 - 7j <J+1 6+1 -6-1 J 

, If fi+l J ~ 

2^ 6+1 -6-1 


(44) 


and the boundary condition is used for the y-component given by 


^Vij- 1/2 



-C ai /3 6 


(45) 
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For a full-span configuration, similar changes are required. For half-node points ( i , J+V 2 , k) that lie 
on the left side surface of the computational box, 


* _]_ ( Vj+l Vj-l * ^t+lj _ Vj+l Vj r 1J-1 

x tj+i/2 2 \ Vj - Tjj—i ^ & +1 - Ci-1 rjj-Vj- i y_1 6+1-ft-l / 

+ l f ^±Li ~ 

2 ^ ft+1-6-1 


(46) 


^Vij+1/2 ^ 



(47) 


At Half-Node Points (i, j, fc±V 2 ) 

At half-node points (z, j, k—V 2 ) only the z-component of the disturbance velocity is required. It 
is defined simply by 


K-m, 


C k - Cfc-i 


(48) 


The formula for the z-component at (z, j, fc+V 2 ) is determined from equation (48) by incrementing the 
k index by one. The z-component of the disturbance velocity, however, needs to be defined differently 
to account for the lifting-surface flow-tangency and wake boundary conditions. The lifting surfaces 
are located an equal distance between grid lines so that in the plane directly above the surface the 
<t>z k _ 1/2 f° rmu l a °f equation (48) is replaced by (/+ + ft) n+1 and in the plane below the surface the 

(j>l k+ 1/2 formula similar to that of equation (48) is replaced by (/“ + ft) n+l - For the wake boundary 
condition the solution procedure is modified by requiring that the disturbance velocity in the vertical 
direction be continuous across the wake (eq. (10b)). This condition is imposed by defining 


*1 - (*U 1 - r j 

Zk ~ l/ 2 Cjt - Cfc-i 


(49) 


at the grid points in the plane above the wake and by defining 

(*Ui- r *)-^ 

^ fc+1 ' 2 Ck+I-Ck 


(50) 


at grid points in the plane below the wake. Furthermore, the z-component of the disturbance 
velocity must be defined differently to account for the flow-tangency boundary condition (eq. (12b)) 
for bodies. The top surface of the computational box that is used to model a body, for example, is 
located equidistantly between grid planes in the vertical direction (i.e., between planes k and k — 1). 
Hence, for grid points in the plane k that is above the top surface, 


j* 

%- 1/2 


- c “(kX- c ^ 


(51) 


Similarly, the bottom surface of the computational box is located equidistantly between grid planes 
in the vertical direction (i.e., between planes k and k+ 1). Hence, for grid points in the plane k that 
is below the bottom surface, 


<t>\ 


z k+ 1/2 


_ r ,N, 

~ - c “ l K 


C ai ab 


(52) 
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Difference Equations for the Grid Metrics 

The finite-difference equations for the grid metrics are presented in this section. These equations 
were derived to be consistent with the differencing of the disturbance potentials, such that if the 
AF algorithm were a full-potential solver rather than a TSD solver, a uniform flow would be an 
exact result of the finite-difference formulas for the velocity components (presented in the previous 
section). The grid metrics are required in numerous places throughout the AF solution procedure 
at grid points (i, j, k) and at the half-node points (iiV 2 , j , k) and (i, jiV 2 , k). 


At Grid Points (i, j , k) 

At grid points (i, j, k), the metrics are defined to be 


t- 


£i+l £i-l 




x i+l,j x i—l,j 
&+1 — &- 1 x jJ ± 1 ~ 


t,] x i+l,j - x i—\,j Vj + 1 “ Vj-1 
However, at the upstream ( i = 1) and downstream (i = NXT) boundaries of the grid, 

. = 1 
Syij= 0 

Along the symmetry plane (j = 1), 

- _ £i+l ~ £i-l x i,j+'t ~ x i,j 

Vl ' j ~ x i+l,j - x i-l,j Vj + 1 - Vj 

and along the far-spanwise boundary ( j = NYT), 

* _ &+1 ~ £i— 1 x ij ~ x ij-l 

^ Vi ' j ~ X i+ lJ - Xi-lj T] j - Vj-l 


( 53 ) 

( 54 ) 

( 55 ) 

( 56 ) 

( 57 ) 


( 58 ) 


At Half-Node Points (z±V 2 , j , k ) 

At the half-node points (z— V 2 , j, k) the metrics are defined to be 




«— 1 / 2,2 


& &— 1 

x i,j ~ x i—l,j 


& ~ &- 1 x iJ + 1 a «-l J+l l,j— 1 

J/i-iAj Xij _ xj.ij 2 (r?j+i - Vj- 1) 

However, along the symmetry plane (j = 1), 


^t-l/2j 


& ~ &-1 g »j+l ~ X hj 4- Xj-lJ+1 x i-l,j 
x i,j x i—l,j 2 (Vj+l Vj) 


and along the far-spanwise boundary (j = NYT), 


tVi-l/2 j 


fjj ~ &-1 x ij ~ x iJ - 1 + x i—l,j Xj-ij-i 
x i,j ~ x i - 1 j 2 “ ^-l) 


( 59 ) 

( 60 ) 


( 61 ) 


( 62 ) 


13 



At Half-Node Points (z, j±V 2 , k ) 

At half-node points (z, j—V: 2 , &), the metrics are defined to be 

^ 2 (&+l -Jj-i) 

^ — ^z— 1 j "1" ^z+l J— 1 + ^z— 1 J— 1 

£ __ ~ 2 te+l ~ £z-l) x i,j ~ ^zj'-l 

J ^ + x z+l J-l + ^z-lj-l Vj ~ Vj-l 


(63) 

(64) 


Difference Equations for the Left-Hand-Side Operators 

The finite-difference equations for the L £, and L ^ operators are presented in this section, 

including the modifications that are required to impose the symmetry plane, lifting-surface flow- 
tangency, and wake boundary conditions and the flow-tangency conditions for bodies. Also, the 
time derivatives in the AF algorithm have been written for variable time stepping to allow for step- 
size cycling to accelerate convergence to steady state. Because of this, three different time steps are 
required as defined by 


At! = f n+1 - t n 

(65) 

A t 2 = t n - t n ~ l 

(66) 

A t 3 = t n ~ l - f n ~ 2 

(67) 


Operator 

The L £ operator is implemented by considering the £-sweep equation defined by 


k 



B 2 Afi + At 2 a d 
2A At!+At 2 


At 1 At 2 , 

2A ■ 


,±Fi±) 

: dt 'od 


A <t> 


= -R 


where, for the isentropic formulation, 


( 68 ) 


£ 2 

F i = E£ x + 2 F&tfe + 2 GZ y <f>* y +fd l + H k) + HZytfy (69) 

for the nonisentropic formulation, 

*i = Zx(l + 1 )Ml 0 QW (V s - V) (1 - K v V) + 2G£ y (f>y + f- (1 + H<f>* x ) + H£ y (f>* (70) 

and R is the residual, which is treated in a subsequent section. The first derivative in the Lg operator 
is represented by a backward-difference formula, to maintain numerical stability, given by 

|( A? ) = s«^: <71) 

The other derivative in the L £ operator is treated by first considering the flux F{ as being the sum 
of two fluxes, so that 


d_ 

9 Z 



9 pEO 9 

a? Fl at 



, 9 P CD 

+ ae F ‘ 


d_ 

9 Z 



(72) 
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where 


( 73 ) 


£ 2 

if D = (2 G + H)t v 4>l + f (1 + ml) 

sx 

and Fp° is taken to be one of two different formulas depending upon whether entropy effects are 
included. For example, in the original isentropic formulation, 

if 0 = EZ X + 2 FU<t>% (74) 


and in the nonisentropic formulation, 

jf 0 = Ul + 1 )M^QW (V s -V)(l- K v V ) (75) 

Regardless of which definition is selected, Fp® is either centrally differenced at subsonic points or 
backward differenced at supersonic points according to the Engquist-Osher (EO) type-dependent 
mixed-difference operator (ref. 13), and Fp® is always centrally differenced (CD) independent of 
the local speed of the flow. Specifically, the second derivative of the first flux is represented by the 
following EO type-dependent formula: 


d „EO d / A 2 IV, . Weo ~ 

dt Fl d£\ A<l> ) i+1/2 >) 1<+ i/2 Hi+i-ti 

J_ ( n . i\ pEO A ^t — A ^t-l „ pE 

+ (tei-1/2 - Ij Flt-ift £ i-V2 F h 


EO A< ^i- 1 ~ A ^t— 2 
■3/2 £i—i — $i—2 . 


where in the isentropic formulation, 


F h-l/2 ~ E &i-l/2j + 2F ^ x i-l/2J^i-l/2j 


'i- 1/2 


= < 


K.1/* * **) 

($**-1/2 J > 


and the sonic streamwise speed (f>x is defined by 


2 F 


(76) 

(77) 

(78a) 

(78b) 

(79) 


p EO 
F b-l/2 


' ’■ r i-l/2j ' 


e t-l/2 - 1 


(^-1/2.2 ^ 
( V i-l/2,j > 
$e»-1/2j 

H/y 1 + 

^ = 1 


(80) 

(81a) 

(81b) 

(83) 

(84) 
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Recall from equations (16a) and (16b) that 


Q = 


2 + ( 7 — ml 


1 1/2 


( 7 + l)M2 

F s = 


Q 2 -l 


2Q 

The derivative involving the second flux (Fp 0 ) is represented by a central-difference formula 
involving terms evaluated at the half-node points in the ^-direction according to 


2 ( F CD A0 i+1 -A^ pCD A <t>i~ A <t>i- 

0( * S{ ) fi+I - {j _, (l+I _ { . V./. • 

^ = < 2G + A-./y + 0 + "«,-,/«) 

sx 2-l/2,j 


where 


(85) 


( 86 ) 


The derivative involving the second flux at the symmetry plane (normally taken to be j = 1), 
however, needs to be defined differently to impose the symmetry condition. This results from 
requiring that 

(87) 

( 88 ) 


Ci/ 2J - = 0 


which implies that 

for grid points in the symmetry plane. 

The £-sweep equation may then be rewritten for solution in quadradiagonal form as 


where 


bi A 4>i_2 + Ci A^j_ 1 + d( A(f> i + ej A = —R^ 
2 1 


u _ ~ Afi A^2 

b’ — £xij 


2A 


6+1 - 6-1 6- 1 ~ 6-2 


•^-3/2^: 


■EO 
L i— 3/2 


(89) 

(90a) 


Ci - - 


+ 


R 2 At i +■ A^2 


2A 


- r, - XI - At 2 &ij 7 V- 

Afi + At 2 !J 6+i - 6-1 


Ati A t 2 . 2 


<XiJ 6+1 - 6-1 L - 6-1 ( 2q " 1/2 + 6-i-6-2 £i - 3/2 ^- 3 / 2 


2A J '6+i -6-1 

Afi A<2 ^ 2 1 pCD 

2^4 te<J 6+ 1 -6- 1 6-6- 1 li-i/2 


(90b) 


dj — 1 + 


£ 2 Ati + At 2 A , , 2 

221 At 1+ Ai 2 2 ^ 6+1 -6-1 


+ 2 A ^ 6 +1 - 6-1 . 6+1 - 6 £i+1/2 ) 6 - 6-1 ( 2£i-1/2 " *) 

_ 2 r 1 pCD | 1 f cd 

2-4 ,J 6+i - 6-i [6+i - 6 li+1/2 6 - 6-i li-1/2 . 


EO 


+ 


■1/2J 
(90c) 
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Ci = 


At x At 2f 2 1 


(i - ^+ 1 / 2 ) ^ 


(90d) 


At\ A<2, 
2 A 1 


B .. £ ± f, cd 

r ^ + i-&_i& +1 -& li+1/2 


Note that it is more efficient computationally to first calculate fy, c^, and e* and then evaluate the 
main diagonal according to 

di = 1 - bi - Ci - ei (91) 

Also, bi = 0 for a purely subsonic flow, so that further computational efficiency may be obtained by 
using a tridiagonal inversion procedure to calculate A (j) rather than the more costly quadradiagonal 
procedure that is necessary for mixed subsonic-supersonic type flows. 

Modifications to the L ^ operator are required to accommodate the flow-tangency boundary 
conditions for bodies. For example, inside the computational box that is used to model a body 
there is, of course, no flow. Therefore, at grid points that lie inside, 


Acj) =0 


(92) 


which results in the following quadradiagonal coefficients and residual (or right-hand side): 


bi = 0 

Ci = 0 
di = 1 
ei = 0 

Ri = 0 


(93a) 

(93b) 

(93c) 

(93d) 

(93e) 


L v Operator 

The Ljf operator is implemented by considering the 77 - sweep equation defined by 

= = (94) 

where 

F 2 = j-(1 + H(j>*) (95) 

The spatial derivative of the L ^ operator is represented by a central-difference formula involving 
terms evaluated at the half-node points in the 77 -direction according to 




Vj+1 - Vj-1 


Fo. 


ij+1/2 


A^+i ~ 

Vj + 1 - Vj 


- f 2 . . 


i, j-1/2 


A<t>j - A<j>j_ x 
Vj ~ Vj-l 


(96) 


where 


F \j-in { 1 + H< f > x iij _ 1 / 2 ) 

1/2 

The 77 -sweep equation may then be rewritten for solution in tridiagonal form as 

Cj A + dj A(f>j + ej A(f>j +1 = A<j> 


(97) 


(98) 
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where 


c j ~ -txij 


Ati A ti 


-F%<. 


dj — 1 + 


2 A r)j + i - r)j-\ r] j - Vj-l W" 1 / 2 

A^i At 2 2 1 

2^4 7/ i+ i - rjj.i rj j+ i - rjj 2 < J +1 / 2 


(99a) 


Ai^ 
?x ‘*> 2.A 


2 1 

- 7?j-i Vj - Vj-i 2 ‘ J_1/2 


j £>X{ j 


Vj + 1 Vj-l Vj Vj 

Aii Ai 2 2 1 


*2,,. 


2^4 r? j+ i - r?j_i r? J+1 - ry 'J+ 1 / 2 


(99b) 

(99c) 


Note that it is more efficient computationally to first calculate cj and ej, and then evaluate the main 
diagonal according to 

( 100 ) 


d j — 1 Cj c j 


A special L v operator is defined to impose the symmetry condition (j = J; normally taken to be 
J = 1). This is accomplished by requiring that 




- 1/2 

The spatial derivative of the L ^ operator is then given by 


J+l/2 


( 101 ) 


hi ( A? ) - 


2 F 1 - A(pj 

VJ + 1 - VJ 2i,J+1/2 ^J+l - VJ 


( 102 ) 


which results in an upper bidiagonal 77-sweep equation at the symmetry plane defined by 

cj = 0 (103a) 

Ati Ato 2 

dj -1+ Zx ij — — A — 7- —^ F 2 u+ m ( 103b ) 


2A (vj+i ~ VJ? M+1/2 
_ Ah Ah 2 
J iJ U tju-n-ru? w 


(103c) 


Modifications to the L ^ operator, similar to the modifications made to the L £ operator, are also 
required to accommodate the flow-tangency boundary conditions for bodies. For grid points that lie 
inside the computational box that is used to model a body, 


A<j> =0 

(104) 

ficients and right-hand side: 


Cj — 0 

(105a) 

dj — 1 

(105b) 

e j= 0 

(105c) 

£ 

II 

O 

(105d) 
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For grid points in the plane j that is adjacent to the right side of the computational box to which 
the body boundary conditions are applied, 


dy( A<t> )j- 1/2 ° 

which results in the following tridiagonal coefficients: 


Cj = 0 

(107a) 

A<i At 2 2 1 

3 XiJ 2 A T)j + l - 7/j_l T)j+ 1 - 7 lj 2 J +1 /2 

(107b) 

, Ati At 2 2 1 p 

6j XiJ 2A 7?i+i - Vj-1 Vj+1 ~ Vj 2j+1/2 

(107c) 


For full-span configurations, a similar modification is required to treat the left side of the computa- 
tional box. For grid points in the plane j that is adjacent to the left side, 

^ / * I \ ft /-I Aft\ 


M 


dy V /j+l/2 


which results in 


, Ati Af 2 2 1 

XiJ 2A T/j+l - 77j_! T)j - Tjj- 1 2 i-V 2 

(109a) 

At! At 2 2 1 

j XiJ 2A T)j +1 - Tjj~\ T)j - rtj-i 2 i~ l l 2 

(109b) 

e j — 0 

(109c) 


Lq Operator 

The Z/£ operator is implemented by considering the £-sweep equation defined by 


where 


r Ail At 2 d d 


f 3 = t 

S£ 


F 3 — ) A<t> = A0 


In the present implementation of the AF algorithm, there is no grid shearing in the vertical direction. 
Consequently, the metric is independent of the C-direction, and the £ x in the numerator of the L ^ 
operator identically cancels the £ x in the denominator of F3. The spatial derivative in the operator 
is thus represented quite simply as 


d 2 / A ^ _ 2 [ A</>fc + i - A( 

d ( 2 ^ Cfc+l - Cfc-l L Cfc+l ~ Cfc 


A</>fc+ 1 ~ A<h _ A<j>k - Affa-i 

Cik+i - Ok Cfc — Cfc-i 


The C-sweep equation may then be rewritten for solution in tridiagonal form as 

c fc Afa - 1 + dfc A <f> k + e k A <f> k+ i = A <f> 
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where 


Cfc = 


Ati A f 2 


1 


dk — 1 + 


2A Cfc+l “ Cfc-l C k ~ Cfc-l 

Ati A i 2 2 1 


e k = -- 


2A Cfc+l - Cfc-l Cfc+l - Cfc 
Ati Af 2 2 1 


+ 


Ati Ai 2 


1 


2A Cfc+i Cfc— l Cfc Cfc— l 


2A Cfc+l - Cfc-l Cfc+l - Cfc 


(114a) 

(114b) 

(114c) 


Again note that it is more efficient computationally to first calculate c j. and and then evaluate 
the main diagonal according to 

dfc = l-Cfc-e fc (115) 


The operator, however, needs to be defined differently to account for the lifting-surface flow- 
tangency and wake boundary conditions. Further, the operator must be modified in both the time- 
linearization and the subiteration step. First consider the modification in the time-linearization step. 
The lifting surfaces are located equidistantly between grid lines so that in the plane directly above 
the surface the derivative in equation (112) represented by 

A<ftfc - 

Cfc - Cfc-l 


is replaced by 

+ ft )"* 1 ~ (ft + ft)" 

and in the plane below the surface the derivative in equation (112) represented by 

A<ftfc+ 1 - A<ft fc 

Cfc+l - Cfc 


is replaced by 

Since these new terms are known quantities, they are brought to the right-hand side of the C-sweep 
equation to create bidiagonal equations. For example, for grid points in the plane above the surface 
the tridiagonal coefficients become 


Cfc =0 

Ati At 2 2 1 

2 A Cfc+l - Cfc-l Cfc+l - Cfc 


e fc = — 


Ati A t 2 


1 


2 A Cfc+l “ Cfc-l Cfc+l - Cfc 

and the corresponding right-hand side of the C-sweep equation is 

2 


A? - Atl A ‘ 2 


(/t+ft) n+1 -(ft+/i)' 


2 A Cfc+l - Cfc-l 

For grid points in the plane below the surface the tridiagonal coefficients become 


Ati At2 2 1 

2 A Cfc+l Cfc— l Cfc Cfc— l 


(116a) 

(116b) 

(116c) 


(117a) 
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_ Ati A t 2 2 1 

2A Cfc+l - Cfc— 1 C k - Cfc-1 


(117b) 


= 0 


(117c) 


and the corresponding right-hand side of the (-sweep equation is 


A <f> 


A ti A f 2 2 

Cfc+l - Cfc-i 


(/,-+/.)" +1 


(fi+ft) 


Now consider the modification in the subiteration step. Here, similar changes are made to 
replace the derivatives in equation (112) by the difference in downwash at time levels n + 1 and * for 
grid points in the planes above and below the lifting surface. This results in the same tridiagonal 
coefficients as derived in equations (116) and (117). However, the right-hand side of the £-sweep 
equation remains simply A<£ since 


(/z + + /t)*= ( 

t , \ n+1 

fx+ft) 

(118a) 

(/x + ft) = ( 

f ,\n+l 

fx + ft) 

(118b) 


The L ^ operator is also modified to account for the wake boundary condition in a way similar to 
that for the flow-tangency condition. To accomplish this, the wake circulation T is first calculated 
from equation (10c), which is equivalent to 


-pn+l I pn+1 

x i- 1/2 h-l/2 


= 0 


(119) 


This equation is discretized with second-order- accurate finite-difference approximations for the space 
and time derivatives, which results in 


„ , l/2 At 1+ A* 2 I? +1 - 17 

S ^ /o A /■ >- e\ I 


Ati 


pn 

1 i 


pTl— 1 


+ 


1 ' 2 ^ Ati + A t 2 A t x 

2 Ati + Af 2 i7 +1 - rjL, At! 


r n — v 

L i - 1 i i-i 


At\ + At 2 

71 — 1 \ 


At 2 


A t\ + At 2 


Ati 


Ati + A t 2 A t 2 


= 0 


(120) 


The equation is solved for the unknown circulation 


Ip 4 , which yields 


pTl+1 = 




1/24 1 


1 2 Ati + A t 2 1 \ 
+ 2 Ati + At 2 Ati/ 



pn-f 1 
x t-l 




if 2Ati + At 2 r? Ati rf-rf- 1 
2 ( Ati + Al 2 AIj Ali + A l 2 At 2 

2 At! + Ai 2 1?+ l - rf_, _ Ali !?_! - i?-, 1 \ 

Ati 4- At 2 Ati Ati + At 2 At 2 J 


( 121 ) 


Through use of equation (121), the circulation is convected downstream by calculating for all 
grid points downstream of the trailing edge, beginning with the trailing-edge value Tt e = — <t> t e * 

The wake circulation is then incorporated within the solution procedure by requiring that the 
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disturbance velocity in the vertical direction be continuous across the wake (eq. (10b)). This 
condition is imposed by defining 



2 

A0fc+i - A fa 

A<f> k - (a^_i - r n+1 + r*) ' 

Cfc+i — Cfc— i 

+ 

l 

Cfc - Cfc— i 


at grid points in the plane above the wake and by defining 


a 2 , 2 (A^fc+i - r n+1 + r*) - A<f> k 

9C 2 (k+i - C*— i 4+1 - Cfc 


Cfc ~ Cfc-i 


(122a) 


(122b) 


at grid points in the plane below the wake. Since the circulation terms are known quantities, they are 
brought to the right-hand side of the C-sweep equation. The tridiagonal coefficients of the £-sweep 
are unchanged consequently, and the right-hand side becomes 


a? - 

v 2A 


^ / pn+1 _ p*\ 

Cfc+i - Cfc-i Cfc - Cfc-i ' ' 

at the grid points in the plane above the wake and 


(123a) 


A $ + 


At! A t 2 2 1 

2A Cfc+i ~ Cfc— l Cfc+l - Cfc 


(r n+1 - r) 


(123b) 


at grid points in the plane below the wake. 

The convection of entropy is governed by the same type of equation as the convection of 
circulation. Namely, the entropy is determined by 


0 a+l 


° x i- 1/2 + %-l/2 


,n+l _ 


= o 


(124) 


Approximating the derivatives of equation (124) in the same way as the derivatives of circulation 
results in 



1 

~ 1/2,i & - 6-i 


1 2Af! Af 2 1 y 1 [ s£. I 

2 Ati + At 2 At ! ) p-V2j & - Si-i 


l/ _2 Ati + At 2 a? Ati sf-s" -1 
+ 2 ^ At i + Af 2 Ati Ati 4 Af 2 A t 2 


2 Ati + At 2 af-l 1 ~ 5 f-l 
Ati "h At 2 Ati 


Ati g f_i a?,] 1 \ 

Ati + Af 2 A< 2 J 


(125) 


Through use of equation (125), the entropy is. convected downstream by calculating s" +1 for all grid 
points downstream of shocks beginning with the value at the shock determined by the Rankine- 
Hugoniot shock jump relation (eq. (19)). 

Modifications to the Lq operator, similar to the modifications made to the and L v operators, 
are also required to accommodate the flow-tangency boundary conditions for bodies. For grid points 
that lie inside the computational box that is used to model a body, 


Acf> = 0 


(126) 
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since there is no flow. This results in the following tridiagonal coefficients and right-hand side: 


c k = 0 

dk = 1 


e k = 0 


R k = 0 


(127a) 

(127b) 

(127c) 

(127d) 


For grid points in the plane k that is above the top surface of the computational box, 


^( A0) Jk _ 1/2 = O 

which results in the following tridiagonal coefficients: 

Cfc = 0 


d k = l + 


Ati A f 2 


1 


(128) 

(129a) 
(129b) 

e , — _ * (129c) 

2 A Cfc+l - Cfc-l Cfc+l - Ck 

Similarly, for grid points in the plane k that is below the bottom surface of the computational box, 

Ati A t 2 2 1 


2^4 Cfc+l Cfc — l Cfc+l Cfc 

Ati A t 2 2 1 


Cfc = -■ 


dfc = 1 + 


efc = 0 


2A Cfc+l “ Cfc-l Cfc “ Cfc-l 

Ati At 2 2 1 

2A Cfc+l ~ Cfc-l Cfc ~ Cfc-l 


(130a) 

(130b) 

(130c) 


These modifications (eqs. (126) to (130)) are made for both the time-linearization and subiteration 
steps of the AF algorithm since the current implementation can only treat steady boundary 
conditions for bodies. 


Difference Equations for the Residual 

The finite-difference equations for the residual are presented in this section. These equations are 
derived by first rewriting the residual in the following general form: 


where 


„ , Ati At 2 ( dgo dgi dg 2 dg 3 \ 

R= 2 a [at + -^ + ~d^ + ~ac) 

(131) 

90 = ~ 

(132a) 

51 - E<j>* + F (</»*) 2 + Gfy) 2 + ^ (1 + H4>%) <t>* y 

(132b) 

92 = T - (1 + H(f >* x ) 4>*y 

SX 

(132c) 
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(132d) 



However, in the nonisentropic formulation, 

9\ = (7 + 1 )M^Q (VV S - ^F 2 ) + G {4>lf + & (1 + H<j>* x ) (P* y (132e) 

Also, the spatial fluxes g\, <72 > and £3, which are in the 77-, and ("-coordinate directions, respectively, 
have been written in terms of the x, y , and z physical coordinates for convenience and computational 
efficiency (as is made evident in the following subsections). 

dgo/dt Term 

The dgo/dt term is treated by considering separately each of the two derivatives that make up 
this term in a manner that is consistent with the treatment of similar terms in the L £ operator. In 
other words, the time term is expressed as 

tr-- u*9 

where the <f>^ and 4>\ t . derivatives are then approximated by finite-difference formulas that account 
for variable time stepping given by 


2 A t 2 - (Afi + A f 2 ) + Ati C 

A t\ + A^2 Afi A<2 

Ati 2 A (A t_2 + A f 3 ) C + At 2 ff~ 2 ' 

Af 2 A^2 + A^3 At2 At3 


* _ At i + A^2 
* Ui ~ A < 2 


(134) 



2 

2 Aix + A f 2 (< 


+ 

i- 1 

1 

& 

1 

Aix + Ai 2 

Aii 

At, 

(•A? -Cl)- 

(r'-ei)" 


Atx + Af 2 At^ 


,(135) 


Note that the spatial derivative in is backward differenced consistent with the treatment of the 
similar term in the L $ operator. Furthermore, for time stepping involving a constant step size (i.e., 
At = At\ = At ‘2 = A £3), equations (134) and (135) simplify to the more familiar formulas given by 


tit, = 


2 # -5# + 4 

At 2 


tin = 


2 3 (g - c) - 4 1 - c 1 ) 


£i+l ~ &-1 


2 At 


(136) 

(137) 


dgi/d£ Term 

The dgi/d£ term is treated by first considering the flux g\ as being the sum of two fluxes, so that 

d$ d£ + d£ < ld8 ) 
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where 


(139) 


9p° = E<p* + F ((f>x) 2 

for the isentropic formulation, 

g f° = ( 7 + l)MlQ (vV s - l -V‘ 

for the nonisentropic formulation, and 

if D = Gte) 2 + ^( l + HfM 


9 1 


€x 


(140) 


(141) 


In equation (139) or equation (140), 0 is either centrally differenced at subsonic points or backward 

differenced at supersonic points according to the Engquist-Osher (EO) type-dependent mixed- 
difference operator, and is always centrally differenced (CD) independent of the local speed 
of the flow. Specifically, the derivative of the first part of the flux is represented by the following EO 
type-dependent formula: 


(l - ^+1/2 ) sF° 1/2 + fei-l/2 “ i) 9i° l/2 


dgf° _ 2 

6+1 - 1 

- £ i-3/29l? 3 / 2 + ( £ i+l/2 “ 2£i- 1/2 + £ i-3/2) 


9 1 


where, in the isentropic formulation, 

n = E ^- m + F 


£ i— 1/2 “ < 


(^i-l/2,j - 

ifxi- 1/2J > 0*) 


and the sonic reference flux yf is defined by 


g{ = E<p s x + F(4> s x ) 2 = 


-E 2 

4F 


In the nonisentropic formulation, 


^° 1/2 = (7 + i)a4qi 

(' / i-l/2^“ - 

1° 

fc-l/W s i" 

£ i—l/2,j - \ 

(Vi-i/ 2j > V‘ 

l 1 


^i-l/2j 


T / _ 

i-V2J l + i^ ; _ 

and the sonic reference flux is defined by 

_L 1 

-m 2 qv 




.s _ 7 + 1 pi/s 
2 


01 = 


(142) 


(143) 


(144) 


(145) 


(146) 
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The derivative of the second part of the flux is represented by a central-difference formula 
involving the flux evaluated at the half-node points in the ^-direction according to 




( n CD _ CD \ 

li_l \ 9 l i+l/ 2 j Sl i- 1/2 j) 


where 


Ci-fl £i- 

= G K-tA#) + i 1 + H(t> * x i- w) ***- 1 / 2 , j 

sx i-l/2j 


(147) 

(148) 


The dgi®/d£ term, however, is defined differently to impose the symmetry condition (normally 
taken to be j — 1 ). This results from requiring that 


= 0 


dg? D n 

J k = ° 


(149) 

(150) 


which implies that 

for grid points in the symmetry plane. 
dg2/dr} Term 

The dg^jdj] term is represented by a central-difference formula involving the flux g<i evaluated at 
the half-node points in the //-direction according to 


992 _ 

2 (a, a , ^ 

(151) 

dg 

T ) j+ 1 - T)j—\ \ 52 M+l/2 9 \j- 1/2 J 

where 

_ ix- (- 1 + ^^j-1/2) 1/2 


92 i,j- 1/2 

(152) 

The 9g2/dr] term, however, is defined differently to impose the symmetry condition ( j = 

J; normally 

taken to be J = 1 ). Here, 

1/2 _ ~^Vi,J+ 1/2 

(153) 

and 

1/2 = ^i.J+1/2 

(154) 

These two conditions imply that 

52 i,J- 1/2 = ~ 92 i,J+ 1/2 

(155) 


so that finally at the symmetry plane 


992 2 

dr) VJ+I - riJ 92i ' J+ 1/2 


(156) 


dgz/d C Term 

The dgz/dC, term is represented by a central-difference formula involving the flux 53 evaluated at 
the half-node points in the (^-direction according to 


9gz 2 / 

9 C Cfc+l - Cfc-1 V 53 «M,fc+i/2 93 iJ,k- 1/2 


(157) 
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where 


1 


(158) 


1/2 


'id, k— 1/2 




-61 

Yz k-lj2 


The dg^/dC term, however, is defined differently when the lifting-surface flow-tangency and wake 
boundary conditions and the flow-tangency conditions for bodies are imposed. This term is defined 
differently for grid points in the planes directly above and below the lifting surface and its wake and 
directly above and below the computational surfaces used to model bodies, since the ^-component 
of the disturbance velocity (<£*) is modified there, as discussed in detail previously. 


Difference Equations for the Far-Field Boundary Conditions 

The finite-difference equations for the far-field boundary conditions are presented in this section. 
These conditions involve the upstream, downstream, far-spanwise, upper, and lower boundaries of 
the grid. All these boundaries except for the upstream boundary are represented by nonreflecting 
conditions. At the upstream boundary the flow is assumed to be free stream, and consequently an 
undisturbed flow condition is prescribed. 


Upstream Boundary 

The upstream boundary condition is implemented during the £-sweep of the AF solution 
procedure to determine values of the intermediate potential 0 on the extreme upstream plane of 
grid points (i = 1). The boundary condition (eq. (13a)) is applied along the extreme grid plane, 
which leads to the trivial equation 

A<£ = 0 (159) 

The equation may be written for solution in quadradiagonal form (eq. (89)) as 


b{ A<^_2 + °i A0*-i + d{ A<^ + e* A <^ +1 — — R{ 

where 


bi = 0 

(160a) 

Cj = 0 

(160b) 

di = 1 

(160c) 

ei = 0 

(160d) 

Ri = 0 

(160e) 


Downstream Boundary 

The downstream boundary condition is implemented during the £-sweep of the AF solution 
procedure to determine values of the intermediate potential 0 on the extreme downstream plane 
of grid points (i = NXT). The condition (eq. (13b)) is applied midway between the extreme and 
adjacent grid planes according to 


P i~ I/ 2 Mi— 1/2 +& <-l/a Mi— 1/2 - ° 


where 



—B | A-l/2 ) 
^-V 2 yjQ- 1/2 / 


(161) 


(162a) 
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/ d2 \ V 2 

A-1/2 - (*A + c ._ i/2 j 

(162b) 

C,- 1/2 = E + 2F 4 ,i l _ w 

(162c) 


With second-order-accurate central-difference and one-sided-difference approximations used for the 
space and time derivatives, respectively, the difference equation for equation (161) becomes 


Ip, 


+ 


+ 6 


i- 1/2 I 


0» ~ 0/-i 

, £t— 1 



| Ati 

V / 

' Ati + At 2 

-i-eA 

Ati 

A^ ; 

Ati “b A^2 

1 =0 



K At 2 ) 

' ei-cT 

A^2 


(163) 


By linearizing consistently with the £-sweep of the AF procedure, the difference equation may be 
rewritten for solution in quadradiagonal form (eq. (89)) as 


b{ A <^_ 2 + c{ A <f>i_i + d{ A ^ + e* A<£ i+1 — — R{ 


where 


and 


bj = 0 


Ci 


2 Ati + A^2 _ . 1 

2 Ati (At i + Af 2 ) i_1/2 ^ X '-V2j £ _ 


J 2Ati + At 2 „ \ , 1 

di ~ 2 At, (Ati + A t 2 ) Fi ~ 1 ! 2 + 6 _ e ._ 1 

e ? - = 0 


/?,• 


a 


i— 1/2 


2 (Ati + Af 2 ) 


2 Ati + A t 2 
Ati 


(^^r+Ci-ei) 


Ati 

At^ 


( 0 ? - 1 + ei 



+0; 


i-l/2 


(164a) 

(164b) 

(164c) 

(164d) 


(164e) 


After the ("-sweep of the AF solution procedure is completed, the potentials on the downstream 
boundary must be recomputed (since the potentials on the adjacent plane of grid points have changed 
values because of the rj- and (-sweeps) according to 


A fc = 


, (Hi 4" A0i_ 1) 

di 


(165) 


Far-Spanwise Boundary 

The far-spanwise boundary condition is implemented during the 7 ]- sweep of the AF solution 
procedure to determine values of the intermediate potential (f) on the extreme spanwise plane of grid 


28 



points (j = NYT). The condition (eq. (13e)) is applied midway between the extreme and adjacent 
grid planes according to 



M j— 1/2 + Mj-l /2 + W j- 1/2 ° 


where 

^j-l /2 [ 4A+ Cj^/2 ) 

C j - l , 2 = E + 2F<t>* Xi ._ l/2 


(166) 


(167a) 

(167b) 


With second-order-accurate central-difference and one-sided-difference approximations used for the 
space and time derivatives, respectively, the difference equation for equation (166) becomes 


1 f[2 A*i + A A t_i 

4 J -i/ 2 ^ [ Ati + At 2 \ Aii J Ati+At 2 \ A t 2 ) 


2 A t_i + A i 2 ( 4>j-i ~ </>"-! \ _ Ati ( 4>j - 1 ~ fj- 1 \ 

At\ + A^2 y J At\ + A^2 y A^2 J 


+ tyij- 


tj+lj - $-1 ,j + jjjj- 1 ~ 


-1J- 


1/2 


2 (Cm -Ci-l) 


Vj - Vj- 1 


(168) 


By linearizing consistently with the 77 -sweep of the AF procedure, the difference equation may be 
rewritten for solution in tridiagonal form as 


where 


and 


Cj A <f>j_i + dj A <j>j + ej A = —Rj 


2 A t\ + At 2 jj 1 

_ 4 At\ (Aii + At 2 ) J-1 / 2 r)j - Vj-\ 

2 Aii + Ai 2 1 

j ~ 4 Ati (Aii + Ai 2 ) J_1/2 Vj ~ Vj- 1 

=0 


flj-i/2 

J 4 (Ati + A t 2 ) 


2 Aii 4“ Ai 2 

Aii 


($ “ 1 ” ^i-i) 






(169) 


(170a) 

(170b) 

(170c) 


(170d) 


After the (-sweep of the AF solution procedure is completed, the potentials on the far-spanwise 
boundary (J = NYT) must be recomputed (since the potentials on the adjacent plane of grid points 
have changed values because of the (-sweep) according to 

A<j>j = -j- (Rj + c j A 0 j_i) (171) 

u j 
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For full-span modeling, a similar boundary condition is imposed to determine values of the 
intermediate potential cf> on the j = 1 plane of grid points. The condition (eq. (13f)) is applied 
between the extreme ( j — 1 ) and adjacent grid planes according to 

2 D J+ 1/ 2 (^)j+i/2 ~ ^'+V2 (^) j + 1/2 ~ (^) j +1/2 ~~ 0 ( 172 ) 

where 


^ J + l /2 - 



(173a) 


^•+1/2 - # + 2F0*. 


,7+1/2 


(173b) 


With second-order-accurate central-difference and one-sided-difference approximations used for the 
space and time derivatives, respectively, the difference equation for equation (172) becomes 


4 D j+ 1/2 


([ 


2 Ati + A t 2 / 0/ - 


At! + At 2 V Af x 


+ 


2 Ati + A t 2 ( 0j+i ~ 0"+i 

Ati + At 2 ( Ati 


! ) 

) 


Ati V 

Ati + At 2 \ At 2 y 

Ati / ‘ffi+i ~ ‘ffi+i \ 
Ati + At 2 y At 2 y 


^« J + l /2 


^i+lj+l ~ ~ <Pj-lj fy+l ~ <t>j 


2 (£i+l ~ £i-l) 


Vj+1 - Vj 


= 0 


(174) 


By linearizing consistently with the 77 -sweep of the AF procedure, the difference equation may be 
rewritten for solution in tridiagonal form (eq. (169)) as 


cj A<^_i + dj A(f)j -J- Cj A<^*_^ — — Rj 


where 


and 


2 At\ + A^ 2 ^ 1 

J 4 Ati (Ati + Af 2 ) J+1 / 2 77 / 4.1 - rjj 

_ 2Ati+ At 2 1 

J 4 Ati (Ati + At 2 ) i+1/2 + 77 / 4-1 - 77 / 


(175a) 

(175b) 

(175c) 


R= £ 2+111 

3 4 (Ati + A t 2 ) 


2 Ati + At 2 

Ati 


(*; - q + f j+1 - ^ ,) 


^(*?-*r‘+« +1 -*3S) 


7 + 1/2 


(175d) 


After the £-sweep of the AF solution procedure is completed, the potentials on the far-spanwise 
boundary ( j = 1 ) must be recomputed (since the potentials on the adjacent plane of grid points 
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have changed values because of the £-sweep) according to 


Afij — — (Rj + ej 


( 176 ) 


Upper Boundary 

The upper boundary condition is implemented during the £-sweep of the AF solution procedure 
to determine values of the potential </> n+1 on the extreme upper plane of grid points (fc = NZT). The 
condition (eq. (13c)) is applied midway between the extreme and adjacent grid planes according to 


where 


2 D k-l/2 (<A 


2 ~ k ~ l ■!* \^ + ) fc - 1/2 + ( 0 C + ) fc — 1/2 


= 0 


D k- 1/2 = ( 4A + 


B 2 


1/2 


'k- 1/2, 


C k-l/2 =E + F (</>*. 

ti+i* ~ Cl,* 


^X itk - Zxij 


&+1 - fi -1 


(177) 


(178a) 

(178b) 

(178c) 


With second-order-accurate central-difference and one-sided-difference approximations used for the 
space and time derivatives, respectively, the difference equation for equation (177) becomes 


1 f[ 2Ati + A Aii ( fkZ_K ] \ 

4 fc “l/2j[ A«i + Ai 2 ( Aii ) At 1 + At 2 \ A t 2 ) 


2 Aii + A t 2 ( 4> k t\ ~ g-i \ _ Aii / </>fc-l~^-l l N \ 
Aii + Ai 2 ( Ai! ) Aix + Ai 2 ( A i 2 / 


+ 


1 - d 

Cfc - Cfc-i 


= o 


(179) 


By linearizing consistently with the (/-sweep of the AF procedure, the difference equation may be 
rewritten for solution in tridiagonal form as 


where 


c k A4> k -\ + d k A <p k + e k A<f> k +i = -R k 

(180) 

2 Aii + Ai 2 1 

Ck ~ 4 Aii (Aii + A i 2 ) fc - 1 / 2 Cfc - Cfc-1 

(181a) 

, 2Aii + Ai 2 n 1 

k 4 Aii (Aii + Ai 2 j k 1/2 Cfc - Cfc-l 

(181b) 

o 

II 

Ai 

(181c) 
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and 


Rk — 


D 


fc — 1/2 


4 (Aij + At 2 ) 


2A ‘^ A ' 2 




+ 0: 


**-1/2 


(181d) 


Lower Boundary 

The lower boundary condition is implemented during the C-sweep of the AF solution procedure 
to determine values of the potential </> n+1 on the extreme lower plane of grid points ( k = 1). The 
condition (eq. (13d)) is applied midway between the extreme and adjacent grid planes according to 


Id 


(*f +, ) m/2 - (^ +1 L V 2 


= 0 


where 


C k+ 1/2 / 


1/2 


D k+ 1/2 - ^ 

c k+i/2 = E + f (4>*. k + 4>* Xi k+1 ) 


<t>Uk ~ *U,k 


<£* =e x . .. 

M &+ 1-&-1 


(182) 

(183a) 

(183b) 

(183c) 


With second-order-accurate central-difference and one-sided-difference approximations used for the 
space and time derivatives, respectively, the difference equation for equation (182) becomes 


2 D k+l/2' 


+ 


/ 2 Ati + A t 2 ( - 4>k+i \ Ati / ^k+l ~ ^fc+1 

\ At i + At 2 y At i J Ati + At 2 y At 2 

2 At 1 + At 2 _ A^ ( 4% -4F 1 

At\ + A t 2 y Ati / Ati + A f 2 y Af 2 


j.n+1 _ rf,n+ 

‘Pfc+l 

Cfc+i - C k 


— o 


(184) 


By linearizing consistently with the (/sweep of the AF procedure, the difference equation may be 
rewritten for solution in tridiagonal form (eq. (180)) as 


where 


c k + d k A (j) k + e k A<p k+1 = -R k 


c k = 0 


— 


0 

(185a) 

2 Ati + Af 2 i 1 

4 At! (At! + At 2 ) k+l > 2 + Ck+l ~ Cfc 

(185b) 

2 Ati + A< 2 1 

4 Ati (Ati + Af 2 ) k+1 / 2 Cfc+i - Cfc 

(185c) 
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and 


Rk 


Dk+l/2 2 Afl "I" A^2 / i j,* j,n\ 

4(At, + A ^) A” (h^-h + l+h-h) 



(185d) 


Concluding Remarks 

A time-accurate approximate-factorization (AF) algorithm has been described for solution of 
the three-dimensional unsteady transonic small-disturbance equation. The AF algorithm consists 
of a time-linearization procedure coupled with a subiteration technique. The algorithm is the basis 
for the CAP-TSD (Computational Aeroelasticity Program-Transonic Small Disturbance) computer 
code, which was developed for unsteady aerodynamic and aeroelastic analyses of realistic aircraft 
configurations. The paper described details on the governing flow equations and boundary conditions, 
with an emphasis on documenting the finite-difference formulas of the AF algorithm. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
November 7, 1991 
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